##R-code file for ``A Boom with Review'', AJPS, July 2012.
##Send questions to colaresi@msu.edu
##File created August 10, 2011

library(foreign)
setwd("/Users/MichaelColaresi/Desktop/Research/RiddleModels/AJPS_RR_Models")
require(pscl)


data1<-read.dta("dyad_analysis_20100309_Revised.dta",convert.underscore=TRUE)
data.mids<-subset(data1,Outcome.SPII.High!=0 & is.na(disno)==FALSE)
data.mids.ties.YEARS<-subset(data1,is.na(disno)==FALSE)
#take this next line out to analyze dispute-years rather than disputes. this is now set for disputes
data.mids.ties.DISPUTES<-subset(data.mids.ties.YEARS, year==endyear)
#next line should be changed to ".YEARS", if dispute years are being analyzed
data.mids.ties<-data.mids.ties.DISPUTES
data.mids.ties.temp<-data.mids.ties[,c("Win.SPII.High","Outcome.SPII.High","name.high","name.low","SPII.diff","cap.diff","majpow.diff","force.diff","threat.diff","mid.log.diff","counter","firstmove.diff","US.inv","prev.disputes", "prev.losses")]
data.mids.ties.nomiss<-na.omit(data.mids.ties.temp)
data.mids.ties.nomiss$Outcome.SPII.High<-data.mids.ties.nomiss$Outcome.SPII.High+2


#calculate number of groups for Dems

J1<-max(as.numeric(as.factor(data.mids.ties.nomiss$name.high)))
J2<-max(as.numeric(as.factor(data.mids.ties.nomiss$name.low)))
name.high1<-as.numeric(as.factor(data.mids.ties.nomiss$name.high))
name.low1<-as.numeric(as.factor(data.mids.ties.nomiss$name.low))

olsConstrRE<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + majpow.diff + mid.log.diff + counter + prev.disputes + prev.losses, data=data.mids.ties.nomiss, x=TRUE, y=TRUE)

forJagsConstrRE <- list(y=olsConstrRE$y, x=apply(olsConstrRE$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrRE$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA, uncomment out and change inits for mu1 and mu2 to run model with no constant
#	mu2=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA, uncomment out and change inits for mu1 and mu2 to run model with no constant
	b0=rep(0,8),
	B0=diag(1E-08,8))

initsConstrRE <- list(list(beta=coef(olsConstrRE)[-1],tau0=1,sigtemp1=1,sigtemp2=1, mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

require(rjags)

fooConstrRE <- jags.model(file="oLogitConstrRE_8VARS.bug", data=forJagsConstrRE, inits=initsConstrRE,n.adapt=5000)
outConstrRE <- coda.samples(fooConstrRE, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=500000, thin=10)
outConstrRE.mcmc<- as.mcmc(outConstrRE)
summary(outConstrRE.mcmc)
save.image(,file="RandomEffectModerate8VAR.RData")
#passes heidel.diag

#p(SPII.coef<0)
1-sum(ifelse(outConstrRE.mcmc[,1]>0,1,0))/length(outConstrRE.mcmc[,1])
#Plot Coefficients
#pdf(file="BOBTCovariates1.pdf")
par(mfrow=c(4,2))
names<-c("Oversight","Capabilities", "Threat", "Major Power", "Other MIDs", "Time Since Prev.", "Prev. Disputes", "Prev. Losses")
yranges<-c(3,3,3,2,2,20,10,2)
#counter is different
for(i in 2:8) {
plot(density(outConstrRE.mcmc[,i]),main=paste(as.character(names[i])), ylim=c(0,yranges[i]), xlab="Parameter")
HPD<-HPDinterval(outConstrRE.mcmc[,i],prob=.80)
lines(HPD, c(-.08,-.08),lwd=4, col="red")
rug(mean(outConstrRE.mcmc[,i]), lwd=6)
abline(v=0)
}
#dev.off()


#Predicted Probabilities
library(boot)
p.SPII<-matrix(NA,nrow=length(outConstrRE.mcmc[,1]),ncol=3)
#changes will be from min to max, p.SPII holds max, p.NOSPII holds min
SPII.change<- max(data.mids.ties.nomiss$SPII.diff) - mean(data.mids.ties.nomiss$SPII.diff)
cap.set<- .5 + sd(data.mids.ties.nomiss$cap.diff)
threat.set<- -.5*sd(data.mids.ties.nomiss$threat.diff)
major.set<- 0-mean(data.mids.ties.nomiss$majpow.diff)
mid.log.set<-log(round(exp(mean(data.mids.ties.nomiss$mid.log.diff))))
counter.set<- -mean(data.mids.ties.nomiss$counter)
prev.m.set<- 0
prev.l.set<- 0-mean(data.mids.ties.nomiss$prev.losses)
#define constant that is equal to mu1-mu2
constant<-outConstrRE.mcmc[,9]-outConstrRE.mcmc[,10]
xB<- constant + outConstrRE.mcmc[,1]*SPII.change + outConstrRE.mcmc[,2]*cap.set +  outConstrRE.mcmc[,3]*threat.set + outConstrRE.mcmc[,4]*major.set + outConstrRE.mcmc[,5]*mid.log.set + outConstrRE.mcmc[,6]*counter.set + outConstrRE.mcmc[,7]*prev.m.set + outConstrRE.mcmc[,8]*prev.l.set

#cap.diff=.4, mid.log.diff=1.4 (which is ln(3))
for (i in 1:length(outConstrRE.mcmc[,1])) {
	p.SPII[i,1] <- inv.logit(outConstrRE.mcmc[i,13] - xB[i])
	p.SPII[i,2] <- inv.logit(outConstrRE.mcmc[i,14] - xB[i]) - inv.logit(outConstrRE.mcmc[i,13] - xB[i])
	p.SPII[i,3] <- 1 - inv.logit(outConstrRE.mcmc[i,14] - xB[i])
}

p.NOSPII<-matrix(NA,nrow=length(outConstrRE.mcmc[,1]),ncol=3)
#data are mean deviated, so this should be too
SPII.change<- -mean(data.mids.ties.nomiss$SPII.diff)
xB<-constant + outConstrRE.mcmc[,1]*SPII.change + outConstrRE.mcmc[,2]*cap.set +  outConstrRE.mcmc[,3]*threat.set + outConstrRE.mcmc[,4]*major.set + outConstrRE.mcmc[,5]*mid.log.set + outConstrRE.mcmc[,6]*counter.set + outConstrRE.mcmc[,7]*prev.m.set + outConstrRE.mcmc[,8]*prev.l.set
for (i in 1:length(outConstrRE.mcmc[,1])) {
	p.NOSPII[i,1] <- inv.logit(outConstrRE.mcmc[i,13] - xB[i])
	p.NOSPII[i,2] <- inv.logit(outConstrRE.mcmc[i,14] - xB[i]) - inv.logit(outConstrRE.mcmc[i,13] - xB[i])
	p.NOSPII[i,3] <- 1 - inv.logit(outConstrRE.mcmc[i,14] - xB[i])
}

#plot, pr(win|SPII==.88)-pr(win|SPII==.01)
#pdf(file="BOBTRepublican.pdf")
par(mfrow=c(3,1))
plot(density(p.SPII[,3]),main="Strong Oversight", xlim=c(0,.2))
abline(h=0)
#HPD<-HPDinterval(as.mcmc(p.SPII[,3]),prob=.90)
CCI<-quantile(p.SPII[,3], probs=c(.05,.96))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]), lwd=6)

plot(density(p.NOSPII[,3]),main="Weak Oversight", xlim=c(0,.2))
abline(h=0)
#HPD<-HPDinterval(as.mcmc(p.NOSPII[,3]),prob=.90)
CCI<-quantile(p.NOSPII[,3], probs=c(.05,.96))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.NOSPII[,3]), lwd=6)

plot(density(p.SPII[,3]-p.NOSPII[,3]), main="Difference (Strong - Weak)", xlim=c(-.1,.2))
abline(h=0)
abline(v=0)
#HPD<-HPDinterval(as.mcmc(p.SPII[,3]-p.NOSPII[,3]),prob=.90)
CCI<-quantile(p.SPII[,3]-p.NOSPII[,3],probs=c(.05,.96))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]-p.NOSPII[,3]), lwd=6)
#dev.off()

#plot random effects
#pdf(file="RandomEffects.pdf")
par(mfrow=c(2,1))
plot(density(outConstrRE.mcmc[,9]), main=expression(paste("Posterior Distribution of ",epsilon[1])), xlim=c(0,3))
abline(v=0)
abline(h=0)
plot(density(outConstrRE.mcmc[,10]), main=expression(paste("Posterior Distribution of ",epsilon[2])), xlim=c(0,3))
abline(v=0)
abline(h=0)
#dev.off()

#Estimate Small model with random effects, constrained.
olsConstrSmallRE<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff, data=data.mids.ties.nomiss, x=TRUE, y=TRUE)
deviate.ols.ConstrSmallRE<- olsConstrSmallRE$x[,-1]-mean(olsConstrSmallRE$x[,-1])
forJagsConstrSmallRE <- list(y=olsConstrSmallRE$y, x=deviate.ols.ConstrSmallRE, 
	N=length(olsConstrSmallRE$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA 
#	mu2=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA
	b0=rep(0,1),
	B0=1E-08)

initsConstrSmallRE <- list(list(beta=coef(olsConstrSmallRE)[-1],tau0=1, sigtemp1=1,sigtemp2=1, mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

#require(rjags)

fooConstrSmallRE <- jags.model(file="oLogitConstrSmallRE.bug", data=forJagsConstrSmallRE, inits=initsConstrSmallRE, n.adapt=5000)
outConstrSmallRE <- coda.samples(fooConstrSmallRE, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=500000, thin=10)
outConstrSmallRE.mcmc<- as.mcmc(outConstrSmallRE)
summary(outConstrSmallRE.mcmc)

#p(SPII.coef<0)
#small model
1-sum(ifelse(outConstrSmallRE.mcmc[,1]>0,1,0))/length(outConstrSmallRE.mcmc[,1])


#Calculate and plot predicted probabilities for small model
p.SPII<-matrix(NA,nrow=length(outConstrSmallRE.mcmc[,1]),ncol=3)
#changes will be from min to max, p.SPII holds max, p.NOSPII holds min
SPII.change<- max(data.mids.ties.nomiss$SPII.diff) - mean(data.mids.ties.nomiss$SPII.diff)
constant<- outConstrSmallRE.mcmc[,2]-outConstrSmallRE.mcmc[,3]
xB<-constant + outConstrSmallRE.mcmc[,1]*SPII.change
#cap.diff=.4, mid.log.diff=1.4 (which is ln(3))
for (i in 1:length(outConstrSmallRE.mcmc[,1])) {
	p.SPII[i,1] <- inv.logit(outConstrSmallRE.mcmc[i,6] - xB[i])
	p.SPII[i,2] <- inv.logit(outConstrSmallRE.mcmc[i,7] - xB[i]) - inv.logit(outConstrSmallRE.mcmc[i,6] - xB[i])
	p.SPII[i,3] <- 1 - inv.logit(outConstrSmallRE.mcmc[i,7] - xB[i])
}

p.NOSPII<-matrix(NA,nrow=length(outConstrSmallRE.mcmc[,1]),ncol=3)
#data are mean deviated, so this should be too
SPII.change<- -mean(data.mids.ties.nomiss$SPII.diff)
xB<- constant + outConstrSmallRE.mcmc[,1]*SPII.change
for (i in 1:length(outConstrSmallRE.mcmc[,1])) {
	p.NOSPII[i,1] <- inv.logit(outConstrSmallRE.mcmc[i,6] - xB[i])
	p.NOSPII[i,2] <- inv.logit(outConstrSmallRE.mcmc[i,7] - xB[i]) - inv.logit(outConstrSmallRE.mcmc[i,6] - xB[i])
	p.NOSPII[i,3] <- 1 - inv.logit(outConstrSmallRE.mcmc[i,7] - xB[i])
}

#pdf(file="BOBTRepublicanSMALL.pdf")
par(mfrow=c(3,1))
plot(density(p.SPII[,3]),main="Strong Oversight", xlim=c(0,.2))
abline(h=0)
CCI<-quantile(p.SPII[,3], probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]), lwd=6)

plot(density(p.NOSPII[,3]),main="Weak Oversight", xlim=c(0,.2))
abline(h=0)
CCI<-quantile(p.NOSPII[,3],probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.NOSPII[,3]), lwd=6)

plot(density(p.SPII[,3]-p.NOSPII[,3]), main="Difference (Strong - Weak)", xlim=c(-.1,.2))
abline(h=0)
abline(v=0)
CCI<-quantile(p.SPII[,3]-p.NOSPII[,3],probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]-p.NOSPII[,3]), lwd=6)
#dev.off()





#######################################################
###Estimate Big model with random effects, constrained
olsConstrBigRE<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + majpow.diff + mid.log.diff + counter + force.diff + firstmove.diff + prev.disputes  + prev.losses, data=data.mids.ties.nomiss, x=TRUE, y=TRUE)
forJagsConstrBigRE <- list(y=olsConstrBigRE$y, x=apply(olsConstrBigRE$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrBigRE$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA 
#	mu2=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA
	b0=rep(0,10),
	B0=diag(1E-08,10))

initsConstrBigRE <- list(list(beta=coef(olsConstrBigRE)[-1], tau0=1, sigtemp1=1,sigtemp2=1,  mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

#require(rjags)

fooConstrBigRE <- jags.model(file="oLogitConstrRE_10VARS.bug", data=forJagsConstrBigRE, inits=initsConstrBigRE, n.adapt=5000)
outConstrBigRE <- coda.samples(fooConstrBigRE, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=750000, thin=10)
outConstrBigRE.mcmc<- as.mcmc(outConstrBigRE)
summary(outConstrBigRE.mcmc)

#p(SPII.coef<0)
#Big Model
1-sum(ifelse(outConstrBigRE.mcmc[,1]>0,1,0))/length(outConstrBigRE.mcmc[,1])


#pred. probs CHANGE THESE TO MATCH MED LARGE MODEL
p.SPII<-matrix(NA,nrow=length(outConstrBigRE.mcmc[,1]),ncol=3)
#changes will be from min to max, p.SPII holds max, p.NOSPII holds min
SPII.change<- max(data.mids.ties.nomiss$SPII.diff) - mean(data.mids.ties.nomiss$SPII.diff)
cap.set<- 0
threat.set<- 0
major.set<- 0-mean(data.mids.ties.nomiss$majpow.diff)
mid.log.set<-log(round(exp(mean(data.mids.ties.nomiss$mid.log.diff))))
counter.set<- 0
force.set<-1-mean(data.mids.ties.nomiss$force.diff)
first.set<-0-mean(data.mids.ties.nomiss$force.diff)
prev.m.set<- 0
prev.l.set<- 0
constant<-outConstrBigRE.mcmc[,11]-outConstrBigRE.mcmc[,12]
xB<-constant + outConstrBigRE.mcmc[,1]*SPII.change + outConstrBigRE.mcmc[,2]*cap.set +  outConstrBigRE.mcmc[,3]*threat.set + outConstrBigRE.mcmc[,4]*major.set + outConstrBigRE.mcmc[,5]*mid.log.set + outConstrBigRE.mcmc[,6]*counter.set + outConstrBigRE.mcmc[,7]*force.set + outConstrBigRE.mcmc[,8]*first.set + outConstrBigRE.mcmc[,9]*prev.m.set + outConstrBigRE.mcmc[,10]*prev.l.set

#cap.diff=.4, mid.log.diff=1.4 (which is ln(3))
for (i in 1:length(outConstrBigRE.mcmc[,1])) {
	p.SPII[i,1] <- inv.logit(outConstrBigRE.mcmc[i,15] - xB[i])
	p.SPII[i,2] <- inv.logit(outConstrBigRE.mcmc[i,16] - xB[i]) - inv.logit(outConstrBigRE.mcmc[i,15] - xB[i])
	p.SPII[i,3] <- 1 - inv.logit(outConstrBigRE.mcmc[i,16] - xB[i])
}

p.NOSPII<-matrix(NA,nrow=length(outConstrBigRE.mcmc[,1]),ncol=3)
#data are mean deviated, so this should be too
SPII.change<- -mean(data.mids.ties.nomiss$SPII.diff)
xB<-constant + outConstrBigRE.mcmc[,1]*SPII.change + outConstrBigRE.mcmc[,2]*cap.set +  outConstrBigRE.mcmc[,3]*threat.set + outConstrBigRE.mcmc[,4]*major.set + outConstrBigRE.mcmc[,5]*mid.log.set + outConstrBigRE.mcmc[,6]*counter.set + outConstrBigRE.mcmc[,7]*force.set + outConstrBigRE.mcmc[,8]*first.set + outConstrBigRE.mcmc[,9]*prev.m.set + outConstrBigRE.mcmc[,10]*prev.l.set
for (i in 1:length(outConstrBigRE.mcmc[,1])) {
	p.NOSPII[i,1] <- inv.logit(outConstrBigRE.mcmc[i,15] - xB[i])
	p.NOSPII[i,2] <- inv.logit(outConstrBigRE.mcmc[i,16] - xB[i]) - inv.logit(outConstrBigRE.mcmc[i,15] - xB[i])
	p.NOSPII[i,3] <- 1 - inv.logit(outConstrBigRE.mcmc[i,16] - xB[i])
}

#plot, pr(win|SPII==.88)-pr(win|SPII==.01)
#pdf(file="BOBTRepublicanBIG.pdf")
par(mfrow=c(3,1))
plot(density(p.SPII[,3]),main="Strong Oversight", xlim=c(0,1))
abline(h=0)
CCI<-quantile(p.SPII[,3], probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]), lwd=6)

plot(density(p.NOSPII[,3]),main="Weak Oversight", xlim=c(0,1))
abline(h=0)
CCI<-quantile(p.NOSPII[,3],probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.NOSPII[,3]), lwd=6)

plot(density(p.SPII[,3]-p.NOSPII[,3]), main="Difference (Strong - Weak)", xlim=c(-.1,1))
abline(h=0)
abline(v=0)
CCI<-quantile(p.SPII[,3]-p.NOSPII[,3],probs=c(.065,.935))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]-p.NOSPII[,3]), lwd=6)
#dev.off()


#Plot Control Coefficients from BIG MODEL
#pdf(file="BOBTCovariatesBig.pdf")
par(mfrow=c(3,3))
names<-c("Oversight","Capabilities", "Threat", "Major Power", "Other MIDs", "Time Since Last MID", "Force", "First Move", "Prev. MIDs", " Prev. Losses")
yranges<-c(3,3,3,2,3,13,2,3,13,2)
#counter is different
for(i in 2:10) {
plot(density(outConstrBigRE.mcmc[,i]),main=paste(as.character(names[i])), ylim=c(0,yranges[i]), xlab="Parameter")
HPD<-HPDinterval(outConstrBigRE.mcmc[,i],prob=.80)
lines(HPD, c(-.02,-.02),lwd=4, col="red")
rug(mean(outConstrBigRE.mcmc[,i]), lwd=6)
abline(v=0)
}
#dev.off()


#######################################################
###Estimate Medium Large model with random effects, constrained
olsConstrMedLargeRE<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + mid.log.diff + counter + force.diff + firstmove.diff, data=data.mids.ties.nomiss, x=TRUE, y=TRUE)
forJagsConstrMedLargeRE <- list(y=olsConstrMedLargeRE$y, x=apply(olsConstrMedLargeRE$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrMedLargeRE$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA 
#	mu2=0,   #COMMENTED OUT BECAUSE THEY ARE PARAMTERS NOT DATA 
	b0=rep(0,7),
	B0=diag(1E-08,7))

initsConstrMedLargeRE <- list(list(beta=coef(olsConstrMedLargeRE)[-1], tau0=1, sigtemp1=1,sigtemp2=1, mu1=0,mu2=0,  .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

#require(rjags)

fooConstrMedLargeRE <- jags.model(file="oLogitConstrRE_7VARS.bug", data=forJagsConstrMedLargeRE, inits=initsConstrMedLargeRE, n.adapt=5000)
outConstrMedLargeRE <- coda.samples(fooConstrMedLargeRE, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=750000, thin=10)
outConstrMedLargeRE.mcmc<- as.mcmc(outConstrMedLargeRE)
summary(outConstrMedLargeRE.mcmc)

#p(SPII.coef<0)
#Med. Large Model
1-sum(ifelse(outConstrMedLargeRE.mcmc[,1]>0,1,0))/length(outConstrMedLargeRE.mcmc[,1])
#Plot Coefficients, COntrols from Med. Large model
pdf(file="BOBTCovariatesMedLarge.pdf")
par(mfrow=c(3,2))
names<-c("Oversight","Capabilities", "Threat", "Other MIDs", "Time Since Prev.", "Force", "First Mover")
yranges<-c(3,3,3,2,15,2,2)
#counter is different
for(i in 2:7) {
plot(density(outConstrMedLargeRE.mcmc[,i]),main=paste(as.character(names[i])), ylim=c(0,yranges[i]), xlab="Parameter")
HPD<-HPDinterval(outConstrMedLargeRE.mcmc[,i],prob=.80)
lines(HPD, c(-.03,-.03),lwd=4, col="red")
rug(mean(outConstrMedLargeRE.mcmc[,i]), lwd=6)
abline(v=0)
}
dev.off()

#Plot of Random effects from MedLarge Model



#pred. probs CHANGE THESE TO MATCH MED LARGE MODEL
#SPII.diff + cap.diff + threat.diff + mid.log.diff + counter + force.diff + firstmove.diff
p.SPII<-matrix(NA,nrow=length(outConstrMedLargeRE.mcmc[,1]),ncol=3)
#changes will be from min to max, p.SPII holds max, p.NOSPII holds min
SPII.change<- max(data.mids.ties.nomiss$SPII.diff) - mean(data.mids.ties.nomiss$SPII.diff)
cap.set<- 0
threat.set<- 0
mid.log.set<- log(round(exp(mean(data.mids.ties.nomiss$mid.log.diff))))
counter.set<- 0
force.set<-1-mean(data.mids.ties.nomiss$force.diff)
first.set<-0-mean(data.mids.ties.nomiss$force.diff)
constant<-outConstrMedLargeRE.mcmc[,8]-outConstrMedLargeRE.mcmc[,9]
xB<-constant + outConstrMedLargeRE.mcmc[,1]*SPII.change + outConstrMedLargeRE.mcmc[,2]*cap.set +  outConstrMedLargeRE.mcmc[,3]*threat.set +  outConstrMedLargeRE.mcmc[,4]*mid.log.set + outConstrMedLargeRE.mcmc[,5]*counter.set + outConstrMedLargeRE.mcmc[,6]*force.set + outConstrMedLargeRE.mcmc[,7]*first.set

#cap.diff=.4, mid.log.diff=1.4 (which is ln(3))
for (i in 1:length(outConstrMedLargeRE.mcmc[,1])) {
	p.SPII[i,1] <- inv.logit(outConstrMedLargeRE.mcmc[i,12] - xB[i])
	p.SPII[i,2] <- inv.logit(outConstrMedLargeRE.mcmc[i,13] - xB[i]) - inv.logit(outConstrMedLargeRE.mcmc[i,12] - xB[i])
	p.SPII[i,3] <- 1 - inv.logit(outConstrMedLargeRE.mcmc[i,13] - xB[i])
}

p.NOSPII<-matrix(NA,nrow=length(outConstrMedLargeRE.mcmc[,1]),ncol=3)
#data are mean deviated, so this should be too
SPII.change<- -mean(data.mids.ties.nomiss$SPII.diff)
xB<- constant + outConstrMedLargeRE.mcmc[,1]*SPII.change + outConstrMedLargeRE.mcmc[,2]*cap.set +  outConstrMedLargeRE.mcmc[,3]*threat.set +  outConstrMedLargeRE.mcmc[,4]*mid.log.set + outConstrMedLargeRE.mcmc[,5]*counter.set + outConstrMedLargeRE.mcmc[,6]*force.set + outConstrMedLargeRE.mcmc[,7]*first.set
for (i in 1:length(outConstrMedLargeRE.mcmc[,1])) {
	p.NOSPII[i,1] <- inv.logit(outConstrMedLargeRE.mcmc[i,12] - xB[i])
	p.NOSPII[i,2] <- inv.logit(outConstrMedLargeRE.mcmc[i,13] - xB[i]) - inv.logit(outConstrMedLargeRE.mcmc[i,12] - xB[i])
	p.NOSPII[i,3] <- 1 - inv.logit(outConstrMedLargeRE.mcmc[i,13] - xB[i])
}

#plot, pr(win|SPII==.88)-pr(win|SPII==.01)
#pdf(file="BOBTRepublicanMedLarge.pdf")
par(mfrow=c(3,1))
plot(density(p.SPII[,3]),main="Strong Oversight", xlim=c(0,.7))
abline(h=0)
CCI<-quantile(p.SPII[,3], probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]), lwd=6)

plot(density(p.NOSPII[,3]),main="Weak Oversight", xlim=c(0,.7))
abline(h=0)
CCI<-quantile(p.NOSPII[,3],probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.NOSPII[,3]), lwd=6)

plot(density(p.SPII[,3]-p.NOSPII[,3]), main="Difference (Strong - Weak)", xlim=c(-.05,.6))
abline(h=0)
abline(v=0)
CCI<-quantile(p.SPII[,3]-p.NOSPII[,3],probs=c(.05,.95))
lines(CCI, c(-.08,-.08),lwd=4, col="red")
rug(mean(p.SPII[,3]-p.NOSPII[,3]), lwd=6)
#dev.off()

#plot(density(p.SPII[,3]/p.NOSPII[,3]), main="Relative Risk (Strong/Weak)", xlim=c(0,19))
#abline(h=0)
#abline(v=1)
#CCI<-quantile(p.SPII[,3]/p.NOSPII[,3],probs=c(.05,.95))
#lines(CCI, c(-.005,-.005),lwd=4, col="red")
#rug(mean(p.SPII[,3]/p.NOSPII[,3]), lwd=6)





#Plot Control Coefficients from BIG MODEL
pdf(file="BOBTCovariatesBig.pdf")
par(mfrow=c(3,3))
names<-c("Oversight","Capabilities", "Threat", "Major Power", "Other MIDs", "Time Since Last Disp.", "Force", "First Move", "Prev. MIDs", " Prev. Losses")
yranges<-c(3,3,3,3,3,10,3,3,10,3)
#counter is different
for(i in 2:10) {
plot(density(outConstrBigRE.mcmc[,i]),main=paste(as.character(names[i])), ylim=c(0,yranges[i]), xlab="Parameter")
HPD<-HPDinterval(outConstrBigRE.mcmc[,i],prob=.80)
lines(HPD, c(-.08,-.08),lwd=4, col="red")
rug(mean(outConstrBigRE.mcmc[,i]), lwd=6)
abline(v=0)
}
dev.off()


#plot random effects from MedLarge
pdf(file="RandomEffectsMedlarge.pdf")
par(mfrow=c(2,1))
plot(density(outConstrMedLargeRE.mcmc[,10]), main=expression(paste("Posterior Distribution of ",epsilon[1])), xlim=c(0,3), xlab="Parameter")
abline(v=0)
abline(h=0)
plot(density(outConstrMedLargeRE.mcmc[,11]), main=expression(paste("Posterior Distribution of ",epsilon[2])), xlim=c(0,3), xlab="Parameter")
abline(v=0)
abline(h=0)
dev.off()

####################################################
####Estimate Alternative Big model with random effects, constrained (Change US.inv for prev.disputes)
#olsConstrBig2RE<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + majpow.diff + mid.log.diff + counter + force.diff + firstmove.diff + prev.disputes, data=data.mids.ties.nomiss, x=TRUE, y=TRUE)
#forJagsConstrBig2RE <- list(y=olsConstrBig2RE$y, x=apply(olsConstrBig2RE$x[,-1], 2, function(x)x-mean(x)), 
#	N=length(olsConstrBig2RE$y),
#	J1=J1,
#	J2=J2,
#	namehigh=name.high1,
#	namelow=name.low1,
#	mu1=0,
#	mu2=0,
#	b0=rep(0,9),
#	B0=diag(1E-08,9))

#initsConstrBig2RE <- list(list(beta=coef(olsConstrBig2RE)[-1], tau0=1, sigtemp1=1,sigtemp2=1,  .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

#require(rjags)

#fooConstrBig2RE <- jags.model(file="oLogitConstrBigRE.bug", data=forJagsConstrBig2RE, inits=initsConstrBig2RE, n.adapt=5000)
#outConstrBig2RE <- coda.samples(fooConstrBig2RE, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=500000, thin=10)
#outConstrBig2RE.mcmc<- as.mcmc(outConstrBig2RE)
#summary(outConstrBig2RE.mcmc)

#p(SPII.coef<0)
#Big Model
#1-sum(ifelse(outConstrBig2RE.mcmc[,1]>0,1,0))/length(outConstrBig2RE.mcmc[,1])

#Plot Control Coefficients from BIG MODEL
#pdf(file="BOBTCovariates3.pdf")
#par(mfrow=c(3,1))
#names<-c("Republican","Capabilities", "Major Power", "Threat", "Other MIDs", "counter", "Force", "First Move", "Prev. Disputes.")
#counter is different
#for(i in 7:9) {
#plot(density(outConstrBig2RE.mcmc[,i]),main=paste(as.character(names[i])), ylim=c(0,6))
#HPD<-HPDinterval(outConstrBig2RE.mcmc[,i],prob=.60)
#lines(HPD, c(-.08,-.08),lwd=4, col="red")
#rug(mean(outConstrBig2RE.mcmc[,i]), lwd=6)
#abline(v=0)
#}
#dev.off()



#Plot Post. Dist. of SPII coefficients across models

pdf(file="BOBTmodelsRepublicanNew.pdf")
par(mfrow=c(4,1))

#Small Model
plot(density(outConstrSmallRE.mcmc[,1]), main="Model 1", xlim=c(-1,5),ylim=c(0,.5), xlab="Oversight Coef.")
abline(v=0)
abline(h=0)
HPD<-HPDinterval(outConstrSmallRE.mcmc[,1], prob=.85)
lines(HPD, c(-.008,-.008),lwd=4, col="red")
rug(mean(outConstrSmallRE.mcmc[,1]), lwd=6)

#Moderate Model
plot(density(outConstrRE.mcmc[,1]), main="Model 2" , xlim=c(-1,5),ylim=c(0,.5), xlab="Oversight Coef.")
abline(v=0)
abline(h=0)
HPD<-HPDinterval(outConstrRE.mcmc[,1], prob=.85)
lines(HPD, c(-.008,-.008),lwd=4, col="red")
rug(mean(outConstrRE.mcmc[,1]), lwd=6)

#Medium Large Model
plot(density(outConstrMedLargeRE.mcmc[,1]), main="Model 3" , xlim=c(-1,5),ylim=c(0,.5), xlab="Oversight Coef.")
abline(v=0)
abline(h=0)
HPD<-HPDinterval(outConstrMedLargeRE.mcmc[,1], prob=.85)
lines(HPD, c(-.008,-.008),lwd=4, col="red")
rug(mean(outConstrMedLargeRE.mcmc[,1]), lwd=6)


#Alternate Big Model (Could also use first Big Model, outConstrBigRE (which has US.inv))
plot(density(outConstrBigRE.mcmc[,1]), main="Model 4" , xlim=c(-1,5),ylim=c(0,.5), xlab="Oversight Coef.")
abline(v=0)
abline(h=0)
HPD<-HPDinterval(outConstrBigRE.mcmc[,1], prob=.85)
lines(HPD, c(-.008,-.008),lwd=4, col="red")
rug(mean(outConstrBigRE.mcmc[,1]), lwd=6)
dev.off()





#######################################
#### Robustness Checks for Online######
######################################

######Just Use MIDS With FORCE########
#subset analysis for only mids where at least one state in the dyad used force (use of force of war)
#create new data that include use of force in hihosta and hihostb
data.mids.ties.temp2<-data.mids.ties[,c("Win.SPII.High","Outcome.SPII.High","name.high","name.low","SPII.diff","cap.diff","majpow.diff","force.diff","threat.diff","mid.log.diff","counter","firstmove.diff","US.inv","hihosta","hihostb", "prev.disputes","prev.losses")]
data.mids.ties.nomiss2<-na.omit(data.mids.ties.temp2)
data.mids.ties.nomiss2$Outcome.SPII.High<-data.mids.ties.nomiss2$Outcome.SPII.High+2

some.force<-ifelse(as.numeric(data.mids.ties.nomiss2$hihosta)>3 | as.numeric(data.mids.ties.nomiss2$hihostb)>3,1,0)
data.mids.ties.nomiss.FORCE<-data.mids.ties.nomiss[some.force==1,]
#REcalculate number of groups for Dems for the subset analysis

J1<-max(as.numeric(as.factor(data.mids.ties.nomiss.FORCE$name.high)))
J2<-max(as.numeric(as.factor(data.mids.ties.nomiss.FORCE$name.low)))
name.high1<-as.numeric(as.factor(data.mids.ties.nomiss.FORCE$name.high))
name.low1<-as.numeric(as.factor(data.mids.ties.nomiss.FORCE$name.low))

olsConstrRE.FORCE<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + majpow.diff + mid.log.diff + counter + prev.disputes + prev.losses, data=data.mids.ties.nomiss.FORCE, x=TRUE, y=TRUE)

forJagsConstrRE.FORCE <- list(y=olsConstrRE.FORCE$y, x=apply(olsConstrRE.FORCE$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrRE.FORCE$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,
#	mu2=0,
	b0=rep(0,8),
	B0=diag(1E-08,8))

initsConstrRE.FORCE <- list(list(beta=coef(olsConstrRE.FORCE)[-1],tau0=1,sigtemp1=1,sigtemp2=1,  mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

require(rjags)

fooConstrRE.FORCE <- jags.model(file="oLogitConstrRE_8VARS.bug", data=forJagsConstrRE.FORCE, inits=initsConstrRE.FORCE,n.adapt=5000)
outConstrRE.FORCE <- coda.samples(fooConstrRE.FORCE, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=500000, thin=10)
outConstrRE.FORCE.mcmc<- as.mcmc(outConstrRE.FORCE)
summary(outConstrRE.FORCE.mcmc)
#save.image(,file="RandomEffectModerateSmallBig.RData")

#p(SPII.coef<0)
1-sum(ifelse(outConstrRE.FORCE.mcmc[,1]>0,1,0))/length(outConstrRE.FORCE.mcmc[,1])

HPDinterval(outConstrRE.FORCE.mcmc[,1], prob=.8)

#########Take out US, still some force############
data.mids.ties.nomiss.noUS<-subset(data.mids.ties.nomiss.FORCE,data.mids.ties.nomiss.FORCE[,"US.inv"]==0)

J1<-max(as.numeric(as.factor(data.mids.ties.nomiss.noUS$name.high)))
J2<-max(as.numeric(as.factor(data.mids.ties.nomiss.noUS$name.low)))
name.high1<-as.numeric(as.factor(data.mids.ties.nomiss.noUS$name.high))
name.low1<-as.numeric(as.factor(data.mids.ties.nomiss.noUS$name.low))

olsConstrRE.noUS<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + majpow.diff + mid.log.diff + counter + prev.disputes + prev.losses, data=data.mids.ties.nomiss.noUS, x=TRUE, y=TRUE)

forJagsConstrRE.noUS <- list(y=olsConstrRE.noUS$y, x=apply(olsConstrRE.noUS$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrRE.noUS$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,
#	mu2=0,
	b0=rep(0,8),
	B0=diag(1E-08,8))

initsConstrRE.noUS <- list(list(beta=coef(olsConstrRE.noUS)[-1],tau0=1,sigtemp1=1,sigtemp2=1,  mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

require(rjags)

fooConstrRE.noUS <- jags.model(file="oLogitConstrRE_8VARS.bug", data=forJagsConstrRE.noUS, inits=initsConstrRE.noUS,n.adapt=5000)
outConstrRE.noUS <- coda.samples(fooConstrRE.noUS, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=500000, thin=10)
outConstrRE.noUS.mcmc<- as.mcmc(outConstrRE.noUS)
summary(outConstrRE.noUS.mcmc)

#p(SPII.coef<0)
1-sum(ifelse(outConstrRE.noUS.mcmc[,1]>0,1,0))/length(outConstrRE.noUS.mcmc[,1])
HPDinterval(outConstrRE.noUS.mcmc[,1], prob=.8)
#Central CI, 75%

save.image(,file="RandomEffectModerateSmallBigRobust1.RData")

########Fatalities Analysis for Robustness
library(MCMCpack)
#cutpoints will be 
data.mids.ties.temp<-data.mids.ties
#use only crisis with use of force (not just show)
some.force<-ifelse(as.numeric(data.mids.ties.temp$hihosta)>3 | as.numeric(data.mids.ties.temp$hihostb)>3,1,0)
data.mids.ties.temp<-data.mids.ties.temp[some.force==1,]

#data to use
data.mids.ties.nomiss.FATAL<-na.omit(data.mids.ties.temp[,c("Outcome.SPII.High","year","name.high","name.low","SPII.diff","cap.diff","majpow.diff","force.diff","threat.diff","mid.log.diff","counter","firstmove.diff","US.inv","hihosta","hihostb","fatal.cat","cat.L","cat.H", "prev.disputes","prev.losses")])


#New cutpoints
cutpoints.diff<-c(-751,-501,-101,-26,-2,1,25,100,500,750)

J1<-max(as.numeric(as.factor(data.mids.ties.nomiss.FATAL$name.high)))
J2<-max(as.numeric(as.factor(data.mids.ties.nomiss.FATAL$name.low)))
name.high1<-as.numeric(as.factor(data.mids.ties.nomiss.FATAL$name.high))
name.low1<-as.numeric(as.factor(data.mids.ties.nomiss.FATAL$name.low))

olsConstrRE.FATAL<-lm(as.numeric(as.character(cat.H)) ~ SPII.diff + cap.diff + majpow.diff + US.inv + prev.disputes + prev.losses , data=data.mids.ties.nomiss.FATAL, x=TRUE, y=TRUE)

forJagsConstrRE.FATAL <- list(y=data.mids.ties.nomiss.FATAL$fatal.cat, x=apply(olsConstrRE.FATAL$x[,-1], 2, function(x)x-mean(x)), 
	N=length(ord.fatal.L),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,
#	mu2=0,
	b0=rep(0,6),
	B0=diag(1E-08,6),
	tau=cutpoints.diff)

initsConstrRE.FATAL <- list(list(beta=coef(olsConstrRE.FATAL)[-1],sigma0=1000,sigtemp1=1,sigtemp2=1,  mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))


require(rjags)

fooConstrRE.FATAL <- jags.model(file="IntervalRE_6VARS.bug", data=forJagsConstrRE.FATAL, inits=initsConstrRE.FATAL,n.adapt=5000)
outConstrRE.FATAL <- coda.samples(fooConstrRE.FATAL, variable.names=c("beta","sigma0","sigtemp1","sigtemp2","mu1","mu2"), n.iter=200000, thin=10)
outConstrRE.FATAL.mcmc<- as.mcmc(outConstrRE.FATAL)
summary(outConstrRE.FATAL.mcmc)

#p(SPII.coef<0)
1-sum(ifelse(outConstrRE.FATAL.mcmc[,1]>0,1,0))/length(outConstrRE.FATAL.mcmc[,1])

summary(outConstrRE.FATAL.mcmc[,1]*.8)
HPDinterval(outConstrRE.FATAL.mcmc[,1]*.8, prob=.8)

save.image(,file="RandomEffectModerateSmallBigRobustInterval.RData")

#Need predicted fatalities from interval regression
#Should be able to just take coefficients and plot them--since they are on the right scale.
#multiply by .8
pdf(file="FatalitiesRepublican.pdf")
plot(density(outConstrRE.FATAL.mcmc[,1]*.8), xlim=c(-200, 300), main="From Low Oversight to High Oversight", xlab="Net Enemy Fatalities")
abline(v=0)
abline(h=0)
HPD<-HPDinterval(outConstrRE.FATAL.mcmc[,1]*.8, prob=.9)
lines(HPD, c(-.00008,-.00008),lwd=8, col="red")
rug(mean(outConstrRE.FATAL.mcmc[,1]*.8), lwd=6)
dev.off()


#####Another Robustness Test
#####Control for Polity Score
#DemocracyScores come from ForAnalysis8_2_22_07.dta, in Secrecy Dilemma folder on Macbook Pro, These are Polity scores
DemScores<-read.dta("DemocracyScores.dta", convert.underscore=TRUE)
#Get rid of previous democ.edit variable, since it is confusion
data.mids.ties.pre <- subset(data.mids.ties, select = -democ.edit)
data.mids.ties.withDem <- merge(data.mids.ties.pre, DemScores, by.x=c("statea","year"), by.y=c("ccode","year"), all.x= TRUE)
data.mids.ties.withDem$democ.edit.a <- data.mids.ties.withDem$democ.edit
data.mids.ties.withDem2 <- subset(data.mids.ties.withDem, select= -c(democ, polity, democ.edit))
data.mids.ties.withDem3 <- merge(data.mids.ties.withDem2, DemScores, by.x=c("stateb","year"), by.y=c("ccode","year"), all.x= TRUE)
data.mids.ties.withDem3$democ.edit.b <- data.mids.ties.withDem3$democ.edit
data.mids.ties.withDem4 <- subset(data.mids.ties.withDem3, select= -c(democ, polity, democ.edit))
data.mids.ties.withDem4$democ.edit.a <- ifelse(is.na(data.mids.ties.withDem4$democ.edit.a), 0, data.mids.ties.withDem4$democ.edit.a)
data.mids.ties.withDem4$democ.edit.b <- ifelse(is.na(data.mids.ties.withDem4$democ.edit.b), 0, data.mids.ties.withDem4$democ.edit.b)

#Create democacy scores for HIGH and LOW SPII
data.mids.ties.withDem4$Democ.diff <- NA
data.mids.ties.withDem4$Democ.diff<-ifelse(data.mids.ties.withDem4$SPII.indexa >= data.mids.ties.withDem4$SPII.indexb,  data.mids.ties.withDem4$democ.edit.a - data.mids.ties.withDem4$democ.edit.b, data.mids.ties.withDem4$democ.edit.b -data.mids.ties.withDem4$democ.edit.a) 


#now clean up and recode outcomes
data.mids.ties.withDem.temp<-data.mids.ties.withDem4[,c("Win.SPII.High","Outcome.SPII.High","name.high","name.low","SPII.diff","cap.diff","majpow.diff","force.diff","threat.diff","mid.log.diff","counter","firstmove.diff","US.inv", "Democ.diff", "prev.disputes", "prev.losses")]
data.mids.ties.withDem.nomiss<-na.omit(data.mids.ties.withDem.temp)
data.mids.ties.withDem.nomiss$Outcome.SPII.High<-data.mids.ties.withDem.nomiss$Outcome.SPII.High+2

#Prepare for BOB-T Model
#Calculate number of dems and non-dems
J1<-max(as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.high)))
J2<-max(as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.low)))
name.high1<-as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.high))
name.low1<-as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.low))

olsConstrREwithDEM<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff  + mid.log.diff + counter + force.diff + firstmove.diff + Democ.diff, data=data.mids.ties.withDem.nomiss, x=TRUE, y=TRUE)

forJagsConstrREwithDEM <- list(y=olsConstrREwithDEM$y, x=apply(olsConstrREwithDEM$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrREwithDEM$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
#	mu1=0,
#	mu2=0,
	b0=rep(0,8),
	B0=diag(1E-08,8))

initsConstrREwithDEM <- list(list(beta=coef(olsConstrREwithDEM)[-1],tau0=1,sigtemp1=1,sigtemp2=1,  mu1=0, mu2=0, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

require(rjags)

fooConstrREwithDEM <- jags.model(file="oLogitConstrRE_8VARS.bug", data=forJagsConstrREwithDEM, inits=initsConstrREwithDEM,n.adapt=5000)
outConstrREwithDEM <- coda.samples(fooConstrREwithDEM, variable.names=c("beta","tau","sigtemp1","sigtemp2","mu1","mu2"), n.iter=750000, thin=10)
outConstrREwithDEM.mcmc<- as.mcmc(outConstrREwithDEM)
summary(outConstrREwithDEM.mcmc)

#SPII controlling for democarcy difference
1-sum(ifelse(outConstrREwithDEM.mcmc[,1]>0,1,0))/length(outConstrREwithDEM.mcmc[,1])

HPDinterval(outConstrREwithDEM.mcmc[,1], prob=.8)
quantile(outConstrREwithDEM.mcmc[,1], probs=c(.05,.95))


#Posterior probability that Democ.diff coefficient is greater than zero and 80 percent HPD
1-sum(ifelse(outConstrREwithDEM.mcmc[,8]>0,1,0))/length(outConstrREwithDEM.mcmc[,7])

HPDinterval(outConstrREwithDEM.mcmc[,8], prob=.8)
quantile(outConstrREwithDEM.mcmc[,8], probs=c(.05,.95))


sink("Output20110614.txt")
print("Small Model")
summary(outConstrSmallRE.mcmc)
quantile(outConstrSmallRE.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrSmallRE.mcmc[,1]>0,1,0))/length(outConstrSmallRE.mcmc[,1])

print("Moderate Model")
summary(outConstrRE.mcmc)
quantile(outConstrRE.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrRE.mcmc[,1]>0,1,0))/length(outConstrRE.mcmc[,1])

print("MedLarge Model")
summary(outConstrMedLargeRE.mcmc)
quantile(outConstrMedLargeRE.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrMedLargeRE.mcmc[,1]>0,1,0))/length(outConstrMedLargeRE.mcmc[,1])

print("Big Model")
summary(outConstrBigRE.mcmc)
quantile(outConstrBigRE.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrBigRE.mcmc[,1]>0,1,0))/length(outConstrBigRE.mcmc[,1])

print("Robust Force Only")
summary(outConstrRE.FORCE.mcmc)
quantile(outConstrRE.FORCE.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrRE.FORCE.mcmc[,1]>0,1,0))/length(outConstrRE.FORCE.mcmc[,1])

print("Robust Non-US Only")
summary(outConstrRE.noUS.mcmc)
quantile(outConstrRE.noUS.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrRE.noUS.mcmc[,1]>0,1,0))/length(outConstrRE.noUS.mcmc[,1])


print("Robust Fatalities") outConstrRE.FATAL.mcmc
summary(outConstrRE.FATAL.mcmc)
quantile(outConstrRE.FATAL.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrRE.FATAL.mcmc[,1]>0,1,0))/length(outConstrRE.FATAL.mcmc[,1])

print("Robust with Democracy Control")
summary(outConstrREwithDEM.mcmc)
quantile(outConstrREwithDEM.mcmc[,1], probs=c(.05,.95))
1-sum(ifelse(outConstrREwithDEM.mcmc[,1]>0,1,0))/length(outConstrREwithDEM.mcmc[,1])
#Posterior probability that Democ.diff coefficient is greater than zero and 80 percent HPD
print("Just democray coefficient now, nonsig")
1-sum(ifelse(outConstrREwithDEM.mcmc[,8]>0,1,0))/length(outConstrREwithDEM.mcmc[,7])
quantile(outConstrREwithDEM.mcmc[,8], probs=c(.05,.95))




####TEST-- DELETE AFTER RUN -- FROM HERE

Prepare for BOB-T Model
#Calculate number of dems and non-dems
J1<-max(as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.high)))
J2<-max(as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.low)))
name.high1<-as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.high))
name.low1<-as.numeric(as.factor(data.mids.ties.withDem.nomiss$name.low))

olsConstrREwithDEM<-lm(as.numeric(Outcome.SPII.High) ~ SPII.diff + cap.diff + threat.diff + majpow.diff + mid.log.diff + counter + force.diff + firstmove.diff + Democ.diff, data=data.mids.ties.withDem.nomiss, x=TRUE, y=TRUE)

forJagsConstrREwithDEM <- list(y=olsConstrREwithDEM$y, x=apply(olsConstrREwithDEM$x[,-1], 2, function(x)x-mean(x)), 
	N=length(olsConstrREwithDEM$y),
	J1=J1,
	J2=J2,
	namehigh=name.high1,
	namelow=name.low1,
	mu1=0,
	mu2=0,
	b0=rep(0,9),
	B0=diag(1E-08,9))

initsConstrREwithDEM <- list(list(beta=coef(olsConstrREwithDEM)[-1],tau0=1,sigtemp1=1,sigtemp2=1, .RNG.seed=1234,.RNG.name="base::Mersenne-Twister"))

require(rjags)

fooTEST <- jags.model(file="TEST.bug", data=forJagsConstrREwithDEM, inits=initsConstrREwithDEM,n.adapt=2000)
outTEST <- coda.samples(fooTEST, variable.names=c("beta","tau","sigtemp1","sigtemp2"), n.iter=10000, thin=1)
outTEST.mcmc<- as.mcmc(outTEST)
summary(outTEST.mcmc)





###DELETE TO HERE




#Make a graph of estimated abilities 
##THIS IS UNFINISHED MAKES SURE IT IS RIGHT!!!!!!!!!!
ability.dist<- matrix(NA, nrow=length(data.mids.ties.nomiss$name.high),ncol=length(outConstrRE.mcmc))

	
data.mids.ties.nomiss$Abilities.low[i]<-quantile(exp(ability.dist[i,]), probs=c(.05))
data.mids.ties.nomiss$Abilities.mean[i]<-mean(exp(ability.dist[i,]))
data.mids.ties.nomiss$Abilities.high[i]<-quantile(exp(ability.dist[i,]), probs=c(.95))
}

#run polr

#Bradley-Terry 2, deleting draws


#####Fatalities D
#B-T without Random Effect
ata analysis Poisson





